package gbench.sandbox.matlib;

import org.junit.jupiter.api.Test;

import java.util.ArrayList;
import java.util.Arrays;
import java.util.Stack;
import java.util.stream.Collectors;
import java.util.stream.Stream;

import gbench.common.matlib.algebra.AlgebraEngine;
import gbench.common.matlib.algebra.num4j.XVec;
import gbench.common.matlib.rdd.BasePair;
import gbench.common.tree.LittleTree.Node;

import static gbench.common.matlib.MatlibCanvas.println;
import static gbench.common.matlib.rdd.BasePair.bp;
import static gbench.common.tree.LittleTree.Term.FT;
import static gbench.common.matlib.algebra.lisp.MyRecord.REC;

/**
 * 牛顿插值
 * 
 * @author Administrator
 *
 */
public class JunitNewton {

    /**
     * 索引节点计算
     * 
     * @param points 数值点
     * @param inode  索引节点
     * @return 索引值
     */
    final Double inode_eval(final XVec<? extends BasePair<? extends Number, ? extends Number>> points,
            final Node<XVec<Integer>> inode) {

        final var v = Stream.of(inode.prop("coef")).map(value -> {
            if (value == null) { // 空值
                final var indices = inode.getValue().toArray(); // 索引节点
                double d = 0d;
                if (inode.isLeaf()) { // 叶子节点
                    d = points.get(indices[0])._2().doubleValue();
                } else { // 非叶子节点
                    final var f0 = inode_eval(points, inode.getChild(0)); // 计算左孩子值
                    final var f1 = inode_eval(points, inode.getChild(1)); // 计算右孩子值
                    final var x0 = points.get(indices[0])._1().doubleValue();
                    final var x1 = points.get(indices[indices.length - 1])._1().doubleValue();
                    d = (f1 - f0) / (x1 - x0);
                } // if
                inode.prop("coef", d); // 设置属性
                return d;
            } else { // 有值
                return (Double) value;
            } // if
        }).findFirst().get();

        return v;
    } //

    /**
     * 牛顿插值
     */
    @Test
    public void interpolate() {
        final XVec<BasePair<Number, Number>> points = XVec.OF(bp(0, 1), bp(2, 2), bp(3, 4), bp(5, 6),
                bp(Math.random(), Math.random())); // 数值节点
        final var indices = XVec.OF(points.size(), i -> i); // 索引
        final var xvex = Stream.iterate(0, i -> i + 1).limit(indices.size())
                .map(i -> indices.slidingS(i + 1, 1, true).collect(XVec.xveclc())).collect(XVec.xveclc());
        final var root = new Node<XVec<Integer>>(xvex.last().get(0));
        final var parents = new ArrayList<Node<XVec<Integer>>>(Arrays.asList(root)); // shangji

        xvex.initialX().rev().forEach(column -> { // 逐列处理
            final XVec<Node<XVec<Integer>>> children = column.fmap(e -> new Node<XVec<Integer>>(e));
            children.slidingX(2, 1, true).foreach((i, cc) -> parents.get(i).addChildren(cc));
            parents.clear();
            parents.addAll(children.list());
        }); // forEach

        root.forEach(inode -> {
            println(" | ".repeat(inode.getLevel()) + inode + "\t" + inode_eval(points, inode));
        });

        println("多项式插值");
        final var stack = new Stack<Node<XVec<Integer>>>();
        stack.add(root);
        final var coefs = new ArrayList<Node<XVec<Integer>>>(); // 多项式系数
        while (!stack.empty()) {
            final var p = stack.pop();
            if (!p.isLeaf()) {
                stack.push(p.getChild(0));
            } // if
            coefs.add(p);
        } // while

        final var newton_interpolate_expression = coefs.stream().map(e -> {
            final var coef = e.dblProp("coef");
            final var dd = e.getValue().delete(e.getValue().size() - 1); // 去除最后一项
            final var _terms = dd.fmap(i -> FT("(x-x$0)", i)).collect(Collectors.joining("*"));
            return _terms.length() > 0 ? FT("$0*$1", coef >= 0 ? coef : FT("neg $0 ", -coef), _terms) : FT("$0", coef);
        }).collect(Collectors.joining("+")); // 生成牛顿插值表达式

        final var engine = new AlgebraEngine();
        println(newton_interpolate_expression); // 显示表达式
        final var node = engine.analyze(newton_interpolate_expression);
        final var xs = points.map(p -> p._1()).collect(XVec.xveclc(true)); // 提取 x 值
        final var _xs = xs.flatMap((_i, _x) -> Stream.of((Object) FT("x$0", (Object) _i), _x))
                .collect(XVec.xveclc(true));

        xs.foreach((i, x) -> { // 基本值的校验
            final var bindings = REC(_xs.append(FT("x", i), x).toArray());
            final var y = node.eval(bindings);
            final var _y = points.get(i)._2().doubleValue();
            println("bindings", bindings, "=====>", "y", y, "_y", _y, "error", y - _y);
        });

        // 计算新的数值
        final var bindings = REC(_xs.append("x", 100).toArray());
        final var v = node.evaluate(bindings);
        println(bindings, "=====>", v);

    }

}
